﻿// mod09ga_prepare.cpp : 此文件包含 "main" 函数。程序执行将在此处开始并结束。
//使用C盘 OSGeo4W的GDAL库快速开发GDAL程序。
//云遥服务器上 探测中心雷达拼图v3.0产品 处理tiff格式程序。 2025/4/16
// Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_23.bin
// 将24层数据转为 500m 间隔的Geotiff格式数据。8km以上的通过插值得到500m间隔层。
// 结果总计32层 500m-16000m
// 层数对应关系
//           原始  结果
// 500-8000  0-15  0-15
// 8500      --    16
// 9000      16    17
// ...
// 15500     --    30
// 16000     23    31
//
// 32层数据量太大了，采用1000米间隔高度层，共15层
//           原始  结果
// 1000      1     0
// 2000      3,5,7,9,11,13     1,2,3,4,5,6
// 8000      15    7
// 9000      16    8
// 16000     23    15

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>
using namespace std;
#include "gdal_priv.h"
#include "bzlib.h"

#define UNICODE
#include "tinydir.h"

struct Header {
    char label[4];
    char version[4];
    int file_bytes = 0;
    short mos_id = 0;
    short coordinate = 0;
    char varname[8];
    char desc[64];
    int blockPos = 0;
    int blockLen = 0;
    int timeZone = 0;
    short yr = 0;
    short mon = 0;
    short day = 0;
    short hr = 0;
    short minu = 0;
    short sec = 0;
    int obsSec = 0;
    short obsDates = 0;
    short genDates = 0;
    int genSec = 0;
    int latmin = 0;
    int lonmin = 0;
    int latmax = 0;
    int lonmax = 0;
    int cx = 0;
    int cy = 0;
    int xsize = 0;
    int ysize = 0;
    int dx = 0;
    int dy = 0;
    short height = 0;
    short compress = 0;
    int numradars = 0;
    int unzipbytes = 0;
    short scale = 0;
    short unused = 0;
    char rgnid[8];
    char units[8];
    char reserved[60];

};

/// 结果会按照文件名字符串从小到大排序。
/// 结果只包括文件名，没有路径信息。
vector<string> scanFiles(string indir, string dateTimeStr) {
    tinydir_dir dir;
    wstring wdirstr = wstring(indir.begin(), indir.end());
    wstring wdtStr = wstring(dateTimeStr.begin(), dateTimeStr.end());
    tinydir_open(&dir, (char*)wdirstr.c_str());
    vector<string> res;
    while (dir.has_next)
    {
        tinydir_file file;
        tinydir_readfile(&dir, &file);
        wstring fname1( &file.name[0], &file.name[strlen(file.name)] );
        if (fname1.find(wdtStr) != string::npos) {
            res.push_back( string(fname1.begin(),fname1.end() ));
        }
        tinydir_next(&dir);
    }
    tinydir_close(&dir);
    std::sort(res.begin(), res.end());
    return res;
}
//筛选出第一个处理时间
vector<string> filterFirstProcess(vector<string> fnameArr) {
    // filename should be like Z_RADA_C_BABJ_20250415013103_P_DOR_ACHN_CAP_20250415_012400_07.bin
    vector<string> res;
    if (fnameArr.size() > 0) {
        string filename0 = fnameArr[0];
        if (filename0.size() == 66 || filename0.size()==64 ) {//64 for cref
            string processDtStr = filename0.substr(14,14);
            for (int i = 0; i < fnameArr.size(); ++i) {
                if (fnameArr[i].find(processDtStr) != string::npos) {
                    res.push_back(fnameArr[i]);
                }
            }
        }
    }
    return res;
}


/// dateTimeStr should be 20250416_010000 format.
/// return 0 is OK.
///  the first part of datetime maybe has two different processing time. use the first process time.
int process32Layers(string inDayDir, string dateTimeStr, string outfile) {

    vector<string> filenameArr = scanFiles(inDayDir, dateTimeStr);
    filenameArr = filterFirstProcess(filenameArr);
    if (filenameArr.size() != 24) {
        std::cout << "file count not 24" << endl;
        return 12;
    }
    CPLStringList opt;
    opt.SetNameValue("COMPRESS", "LZW");
    const int outxsize = 6200;
    const int outysize = 4200;
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* dsout = driver->Create(outfile.c_str(),
        outxsize, outysize, 32, GDT_Int16, opt.List());
    if (dsout == 0) {
        std::cout << "failed to create output " << outfile << endl;
        return 13;
    }
    vector<unsigned char> lastDatabuffer(outxsize * outysize * 2);
    int bandIndex0Based = 0;
    for (int i = 0; i < filenameArr.size() ; ++i) {

        string filenameX = inDayDir+"/"+filenameArr[i] ;
        //sprintf(filename, "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_%02d.bin", i);
        FILE* pf = fopen(filenameX.c_str(), "rb");
        if (pf == 0) {
            std::cout << "failed to open " << filenameX << endl;
            return 14;
        }

        Header header1;
        fread(&header1, sizeof(Header), 1, pf);
        int     bzerror = 0;
        if (header1.xsize != outxsize) {
            std::cout << "wrong header1.xsize " << header1.xsize << endl;
            return 15;
        }
        if (header1.ysize != outysize) {
            std::cout << "wrong header1.ysize " << header1.xsize << endl;
            return 16;
        }
        //the pf is already move to 256 , so no need fseak.
        //fseek(pf, 256, SEEK_SET);
        BZFILE* bf = BZ2_bzReadOpen(&bzerror, pf, 0, 0, NULL, 0);
        if (bzerror != BZ_OK) {
            BZ2_bzReadClose(&bzerror, bf); bf = 0;
            fclose(pf); pf = 0;
            std::cout << "BZ2_bzReadOpen is failed." << endl;
            return 17;
        }
        vector<unsigned char> databuffer(header1.xsize * header1.ysize * 2);//data is short.
        int nread = BZ2_bzRead(&bzerror, bf, databuffer.data(), databuffer.size());
        if (bzerror != BZ_OK && bzerror != BZ_STREAM_END) {
            //BZ_STREAM_END is OK, not error.
            std::cout << "bad  " << bzerror << endl;
        }
        if (bf != 0) BZ2_bzReadClose(&bzerror, bf); bf = 0;
        if (pf != 0) fclose(pf); pf = 0;

        int asize = outxsize * outysize;
        short* sdata = (short*)databuffer.data();
        for (int it = 0; it < asize; ++it) {
            if (sdata[it] < 0) {
                sdata[it] = 0;//use 0 to fill all negative values.
            }
        }
        std::cout << i << " height " << header1.height << " scale " << header1.scale << endl;
        if (i <= 15) {
            dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0, outxsize, outysize, databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
            bandIndex0Based++;
        }
        else {
            //interpol between heights.
            short* dataLow = (short*)lastDatabuffer.data();
            short* dataHigh = (short*)databuffer.data();
            for (int it = 0; it < asize; ++it) {
                if (dataLow[it] > 0 && dataHigh[it] > 0) {
                    dataLow[it] = (dataLow[it] + dataHigh[it]) / 2.0;
                }
                else {
                    dataLow[it] = 0;//the value lower than 0 are all fill with 0.
                }
            }
            dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0,
                outxsize, outysize, lastDatabuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
            bandIndex0Based++;
            dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0,
                outxsize, outysize, databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
            bandIndex0Based++;
        }
        std::cout << "bandIndex0Based " << bandIndex0Based << endl;
        if (i >= 15) {
            lastDatabuffer = databuffer;
        }
    }
    double trans[] = { 73.0, 0.01, 0.00 , 54.2, 0.00, -0.01 };
    dsout->SetGeoTransform(trans);     // trans[6]
    dsout->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
        "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,"
        "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //使用这个不需要注册环境变量，也不需要引用ogr的头文件
    GDALClose(dsout);
    return 0;
}

// 1000      1     0
// 2000      3,5,7,9,11,13     1,2,3,4,5,6
// 8000      15    7
// 9000      16    8
// 16000     23    15
int process16Layers(string inDayDir, string dateTimeStr, string outfile) {

    vector<string> filenameArr = scanFiles(inDayDir, dateTimeStr);
    filenameArr = filterFirstProcess(filenameArr);
    if (filenameArr.size() != 24) {
        std::cout << "file count not 24" << endl;
        return 12;
    }
    CPLStringList opt;
    opt.SetNameValue("COMPRESS", "LZW");
    const int outxsize = 6200;
    const int outysize = 4200;
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* dsout = driver->Create(outfile.c_str(),
        outxsize, outysize, 16, GDT_Int16, opt.List());
    vector<int> usedLayers = { 1,3,5,7,9,11,13,15,16,17,18,19,20,21,22,23 };
    if (dsout == 0) {
        std::cout << "failed to create output " << outfile << endl;
        return 13;
    }
    int bandIndex0Based = 0;
    for (int i = 0; i < filenameArr.size(); ++i) {
        bool isUse = false;
        for (int j = 0; j < usedLayers.size(); ++j) {
            if (usedLayers[j] == i) {
                isUse = true;
                break;
            }
        }
        if (isUse == false)continue;//skip
        string filenameX = inDayDir + "/" + filenameArr[i];
        //sprintf(filename, "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_%02d.bin", i);
        FILE* pf = fopen(filenameX.c_str(), "rb");
        if (pf == 0) {
            std::cout << "failed to open " << filenameX << endl;
            return 14;
        }

        Header header1;
        fread(&header1, sizeof(Header), 1, pf);
        int     bzerror = 0;
        if (header1.xsize != outxsize) {
            std::cout << "wrong header1.xsize " << header1.xsize << endl;
            return 15;
        }
        if (header1.ysize != outysize) {
            std::cout << "wrong header1.ysize " << header1.xsize << endl;
            return 16;
        }
        //the pf is already move to 256 , so no need fseak.
        //fseek(pf, 256, SEEK_SET);
        BZFILE* bf = BZ2_bzReadOpen(&bzerror, pf, 0, 0, NULL, 0);
        if (bzerror != BZ_OK) {
            BZ2_bzReadClose(&bzerror, bf); bf = 0;
            fclose(pf); pf = 0;
            std::cout << "BZ2_bzReadOpen is failed." << endl;
            return 17;
        }
        vector<unsigned char> databuffer(header1.xsize * header1.ysize * 2);//data is short.
        int nread = BZ2_bzRead(&bzerror, bf, databuffer.data(), databuffer.size());
        if (bzerror != BZ_OK && bzerror != BZ_STREAM_END) {
            //BZ_STREAM_END is OK, not error.
            std::cout << "bad  " << bzerror << endl;
        }
        if (bf != 0) BZ2_bzReadClose(&bzerror, bf); bf = 0;
        if (pf != 0) fclose(pf); pf = 0;

        int asize = outxsize * outysize;
        short* sdata = (short*)databuffer.data();
        for (int it = 0; it < asize; ++it) {
            if (sdata[it] < 0) {
                sdata[it] = 0;//use 0 to fill all negative values.
            }
        }
        std::cout << i << " height " << header1.height << " scale " << header1.scale << endl;
        dsout->GetRasterBand(bandIndex0Based + 1)->RasterIO(GF_Write, 0, 0, outxsize, outysize,
            databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);
        bandIndex0Based++;
    }
    double trans[] = { 73.0, 0.01, 0.00 , 54.2, 0.00, -0.01 };
    dsout->SetGeoTransform(trans);     // trans[6]
    dsout->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
        "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,"
        "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //使用这个不需要注册环境变量，也不需要引用ogr的头文件
    GDALClose(dsout);
    return 0;
}

//one layer
int processCref(string inDayDir, string dateTimeStr, string outfile) {

    vector<string> filenameArr = scanFiles(inDayDir, dateTimeStr);
    filenameArr = filterFirstProcess(filenameArr);
    if (filenameArr.size() == 0) {
        std::cout << "file count == 0" << endl;
        return 12;
    }
    CPLStringList opt;
    opt.SetNameValue("COMPRESS", "LZW");
    const int outxsize = 6200;
    const int outysize = 4200;
    GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
    GDALDataset* dsout = driver->Create(outfile.c_str(),
        outxsize, outysize, 1, GDT_Int16, opt.List());
    if (dsout == 0) {
        std::cout << "failed to create output " << outfile << endl;
        return 13;
    }

    string filenameX = inDayDir + "/" + filenameArr[0];
    //sprintf(filename, "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_%02d.bin", i);
    FILE* pf = fopen(filenameX.c_str(), "rb");
    if (pf == 0) {
        std::cout << "failed to open " << filenameX << endl;
        return 14;
    }

    Header header1;
    fread(&header1, sizeof(Header), 1, pf);
    int     bzerror = 0;
    if (header1.xsize != outxsize) {
        std::cout << "wrong header1.xsize " << header1.xsize << endl;
        return 15;
    }
    if (header1.ysize != outysize) {
        std::cout << "wrong header1.ysize " << header1.xsize << endl;
        return 16;
    }
    //the pf is already move to 256 , so no need fseak.
    //fseek(pf, 256, SEEK_SET);
    BZFILE* bf = BZ2_bzReadOpen(&bzerror, pf, 0, 0, NULL, 0);
    if (bzerror != BZ_OK) {
        BZ2_bzReadClose(&bzerror, bf); bf = 0;
        fclose(pf); pf = 0;
        std::cout << "BZ2_bzReadOpen is failed." << endl;
        return 17;
    }
    vector<unsigned char> databuffer(header1.xsize * header1.ysize * 2);//data is short.
    int nread = BZ2_bzRead(&bzerror, bf, databuffer.data(), databuffer.size());
    if (bzerror != BZ_OK && bzerror != BZ_STREAM_END) {
        //BZ_STREAM_END is OK, not error.
        std::cout << "bad  " << bzerror << endl;
    }
    if (bf != 0) BZ2_bzReadClose(&bzerror, bf); bf = 0;
    if (pf != 0) fclose(pf); pf = 0;

    int asize = outxsize * outysize;
    short* sdata = (short*)databuffer.data();
    for (int it = 0; it < asize; ++it) {
        if (sdata[it] < 0) {
            sdata[it] = 0;//use 0 to fill all negative values.
        }
    }
    dsout->GetRasterBand(1)->RasterIO(GF_Write, 0, 0, outxsize, outysize,
        databuffer.data(), outxsize, outysize, GDT_Int16, 0, 0, 0);

    double trans[] = { 73.0, 0.01, 0.00 , 54.2, 0.00, -0.01 };
    dsout->SetGeoTransform(trans);     // trans[6]
    dsout->SetProjection("GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],"
        "AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,"
        "AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]");     //使用这个不需要注册环境变量，也不需要引用ogr的头文件
    GDALClose(dsout);
    return 0;
}



int main(int argc,char* argv[])
{
    std::cout << "v0.1 A program to merge weather radar data into geotiff." << endl;
    std::cout << "v0.2 add cref " << endl;
    std::cout << "usage:yy_radar [mode] /some/dir datetime /some/out/filename.tif" << endl;
    std::cout << "mode: cap/cref , default is cap" << endl;
    std::cout << "datetime should be in format like '20250416_010000' " << endl;
    if (argc != 4 && argc != 5) {
        std::cout << "argc not 4" << endl;
        return 11;
    }
    GDALAllRegister();
    string mode = "cap";
    string inDayDir  ;
    string datetimeStr  ;
    string outfile ;
    if (argc == 4) {
        mode = "cap";
        inDayDir = argv[1];
        datetimeStr = argv[2];
        outfile = argv[3];
    }else if (argc == 5) {
        mode = argv[1];
        inDayDir = argv[2];
        datetimeStr = argv[3];
        outfile = argv[4];
    }

    //int resultCode = process32Layers(inDayDir, datetimeStr, outfile );
    if (mode.compare("cap") == 0) {
        int resultCode = process16Layers(inDayDir, datetimeStr, outfile);
        std::cout << "resultCode " << resultCode << endl;
    }
    else if (mode.compare("cref") == 0) {
        int resultCode = processCref(inDayDir, datetimeStr, outfile);
        std::cout << "resultCode " << resultCode << endl;
    }
    
    std::cout << "done" << endl;
    return 0;
}


/*
* CAP是分层的 后面的数字对应的事高度，8km以下是500m一层，8km以上是1km一层，一共出了24层，对 是插值到某个高度层的反射率，对 0对应500m 1对应1km
0	500
1	1000 0
2	1500
3	2000 1
4	2500
5	3000 2
6	3500
7	4000 3
8	4500
9	5000 4
10	5500
11	6000 5
12	6500
13	7000 6
14	7500 
15	8000 7
16	9000 8
17	10000 9
18	11000 10
19	12000 11
20	13000 12
21	14000 13
22	15000 14
23	16000 15
    vector<string> infiles = {
        "E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_20.bin"
        ,"E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_21.bin"
        ,"E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_22.bin"
        ,"E:/TempAndClean/Z_RADA_C_BABJ_20250416010700_P_DOR_ACHN_CAP_20250416_010000_23.bin"
        ,"E:/TempAndClean/Z_RADA_C_BABJ_20250302001907_P_DOR_ACHN_CAP_20250302_001200_06.bin"
           ,"E:/TempAndClean/Z_RADA_C_BABJ_20250302001907_P_DOR_ACHN_CAP_20250302_001200_11.bin"
    };
for (int i = 0; i < infiles.size() ; ++i) {
        FILE* pf = fopen(infiles[i].c_str() , "rb");
        if (pf == 0) {
            cout << "failed to open " << infiles[i] << endl;
            return 12;
        }
        Header header1;
        fread(&header1, sizeof(Header), 1, pf);

        fclose(pf);
        cout << i << " height " << header1.height << endl;
    }

*/